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Abstract 

The annual modulation in the rate of WIMP recoils observed by the DAMA collaboration at 
high significance is often analyzed in the context of an isothermal Maxwell-Boltzmann velocity 
distribution. While this is the simplest model, there is a need to consider other well motivated 
theories of halo formation. In this paper, we study a different halo model, that of self similar 
infall which is characterized by the presence of a number of cold streams and caustics, not seen in 
simulations. It is shown that the self similar infall model is consistent with the DAMA result both 
in amplitude and in phase, for WIMP masses exceeding ~ 250 GeV at the 99.7% confidence level. 
Adding a small thermal component makes the parameter space near m x = 12 GeV consistent 
with the self similar model. The minimum x 2 per degree of freedom is found to be 0.92(1.03) 
with(without) channeling taken into account, indicating an acceptable fit. For WIMP masses much 
greater than the mass of the target nucleus, the recoil rate depends only on the ratio a p /m x which 
is found to be ~ 0.06 femtobarn/TeV. However as in the case of the isothermal halo, the allowed 
parameter space is inconsistent with the null result obtained by the CDMS and Xenon experiments 
for spin-independent elastic scattering. Future experiments with directional sensitivity and mass 
bounds from accelerator experiments will help to distinguish between different halo models and/or 
constrain the contribution from cold flows. 
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I. INTRODUCTION 



It was shown by Drukier, Freese, and Spergel [1], and by Freese, Frieman, and Gould 
[2] that the motion of the earth about the sun introduces an annual modulation in the flux 
of dark matter particles reaching the earth. The detection of such an annual modulation 
has been claimed by the DAMA/Nal and DAMA/LIBRA experiments 0, E] conducted 
at the Gran Sasso National Laboratory using highly pure Nal(Tl) detectors. The DAMA 
experiment has reported its results for a cumulative time period of 13 annual cycles and a 
total exposure of 1.17 ton- year, claiming a detection of the annual modulation signature at 
> 8cr [3 J . The DAMA claim is strengthened by the fact that only single hit events (expected 
to be triggered by particles with a weak cross section) are annually modulated, the multiple 
hit events show no statistically significant modulation. We refer the reader to (31 H] for 
details regarding the experimental setup and backgrounds. 

The annual modulation seen by the DAMA experiment is commonly analyzed in the 
context of an isothermal Maxwell-Boltzmann velocity distribution implying a WIMP mass 
m x pd 12 GeV or m x ~ 78 GeV for the simple case of elastic spin-independent scattering. 
The derived values of mass and cross section are inconsistent with the null result obtained 
by other dark matter direct detection experiments such as CDMS [5 J and Xenon [6], for 
the simple case of spin-independent elastic scattering. The low mass region may also be 
challenged by observations of the CMB [THS], or by future accelerator experiments. 

The purpose of this paper is to compare the DAMA results to a non-standard halo model, 
namely self similar infall. The self similar infall halo is characterized by a number of discrete 
cold flows and caustics, not seen in numerical simulations. The presence of a cold flow is 
significant since dark matter detection experiments such as DAMA are sensitive to the local 
phase space distribution. The annual modulation effect predicted by the self similar infall 
halo model was studied by Copi and Krauss [10], Green [11], Gelmini and Gondolo [T2] . 
Vergados [15] . and Ling, Sikivie, and Wick [H]. In these papers, it was shown that the self 
similar model predicts qualitatively different results than those predicted by the Maxwellian 
halo. However this does not mean that the maximum recoil rate observed by DAMA on 
May 25 ± 8 days in the 2 — 6 keV ee range [I] is inconsistent with the self similar infall 
model. We show here that for WIMP masses exceeding 250 GeV, the self similar model is 
in agreement with the DAMA observation at the 99.7% level. Nevertheless as in the case 
of the isothermal halo, the allowed parameter space is in contradiction with the exclusion 
limits obtained by other experiments for spin-independent elastic scattering. 

In section II, we derive the recoil rate observed by the DAMA experiment in terms of the 
mass, cross section and velocity distribution. We then briefly discuss the self similar infall 
halo model in Section III. The model is characterized by a series of cold flows, one of which 
is dominant due to the presence of a nearby dark matter caustic. The fractional density 
contributed by the dominant flow is fixed by requiring that the recoil rate be a maximum 
on the observed date of May 25 ± 8 days. We discuss the effect of the dominant flow on 
the annual modulation for Na and I nuclei. We then present our results in Section IV. We 
present the spectrum of expected recoil events for different times of the year, for 4 different 
energy bins. A x 2 analysis is performed to determine the best fit values of mass and cross 
section. We show that the \ 2 P er degree of freedom is close to unity, indicating an acceptable 
fit. We compare the results with that of a Maxwellian halo, and also with the CDMS and 
Xenon bounds. We then check that introducing a small thermal component does not lead 
to qualitatively different predictions. Finally, we verify that the time averaged recoil rate 
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predicted by the self-similar model is consistent with the average plus background reported 
by the DAMA collaboration. In Section V, we summarize our results and discuss various 
ways in which the model may be tested by future experiments. 



II. RECOIL RATE. 




FIG. 1: WIMP-nucleus scattering. 

Consider an elastic collision between a dark matter particle with mass m x and a target 
nucleus of mass m^. The dark matter particle has a velocity v = vx relative to the target 
nucleus. The velocity of the center of momentum is given by v cm = m "^^ N X. The velocity 
of the recoiling nucleus in the CM frame is ^.cm = ~ m m ^ N »- Therefore, in the lab frame 
(where the detector is at rest), the recoil velocity is t>N jlab = "Wn.cm + ^cm = m m +m N (& ~ ""-)• 
The kinetic energy of the recoiling nucleus in the lab frame is: 

Q = = — =^(1 -cos0) , (1) 

2 77T.N 

where 9 is the scattering angle in the CM frame, and = m x m^/ {m x + m^) is the WIMP- 
nucleus reduced mass. The maximum possible recoil energy when the WIMP has a speed v 
relative to the detector is obtained when 9 = n: 

Qmax j (2) 

m N 

and therefore, the minimum velocity the WIMP must have in order to effect a recoil at 
energy Q is found to be 

For example, a 100 GeV WIMP moving at a speed 10~ 3 c and colliding with a w 120 GeV 
Iodine nucleus deposits energy ps 25 keV (1 — cos^). 

The number of recoil events seen by the detector per unit time per unit detector mass 
and per unit recoil energy is 

dR 1 p x da 

— = — / — vf(v)dv. (4) 

dQ m N m x Jv min (Q) dQ 



3 



a is the WIMP-nucleus scattering cross section, p x is the dark matter density at the earth's 
location, and f(v) is the speed distribution of WIMPs relative to the detector. The differ- 
ential cross section is commonly expressed as: 



da (T ,, 2 



dQ Q T 



F (Q) 



2m 2 l v 2 



F\Q). (5) 



F(Q) is called the form factor and contains the momentum dependence of the cross section. 
We assume F(Q) may be described by the form [ToT - fTT] 

F(Q) = e -|M 2 , (6) 

qr 

in units where h and c are set to 1. q = \/2Qm^^ s = 1 fm, R = 1.2 A 1 / 3 fm, r = \/R 2 — 5s 2 , 
and ji is the spherical bessel function. F 2 (Q) w 1 for small Q, and falls off at large Q. We 
will also write <7o in terms of the scattering cross section with a proton or neutron o~ p : 

a = a p A'(^-)\ (7) 

A is the atomic mass number and m^p = m x m p / (m x + m p ) is the WIMP-proton reduced 
mass (we ignore the neutron-proton mass difference). We only consider spin-independent 
elastic scattering here. The recoil rate is 

T is the mean inverse speed 

T=r dv f M (9) 

which is time dependent due to the earth's motion about the sun. It is this term that leads 
to the annual modulation in recoil energy. 

In order to apply Eq ^ to the DAMA experiment, we need to take into account (i) only 
a fraction of the deposited energy is detected and (ii) the target makes use of 2 elements, 
namely Na and I. The fraction of the recoil energy detected by the DAMA experiment is 
called the quenching factor which we label by q.f.(X) where X could be Na or I. In order to 
write Eq ^ in terms of the detected energy Qdet, we make use of the equality of the total 
number of events 

dR . „ dR . — ^ 

— AQ = —— AQ det , (10) 



and therefore, 



dR . . dR 
dQZt ' = dQ 



AQ 



det 



q.f.(A) dQ' 



dR (id 



We can now write down the complete formula for the recoil rate: 

^ ^-(Na) + — ^-r-^(I). (12) 



dQdct Af<j a + A\ dQdet ^Na + Ai dQdet 
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We express our energies in terms of Qdet since this is the quantity measured by the exper- 
iment. We use the unit "keV electron-equivalent (keV ee )" to indicate that Qdet is being 
measured. 

Since the DAMA experiment uses a crystalline detector, it is possible for an ion or recoiling 
nucleus moving parallel to the crystal axes to penetrate deep into the material. Such an ion 
is said to be channeled. Channeled ions transfer their energy primarily to electrons, leading 
to a near unity quenching factor [18], i.e. Q = Qdet- When the channeled fraction is known, 



the effect can be included in the calculation of the total recoil rate Eq (12). 



III. THE SELF SIMILAR INFALL HALO MODEL 



A galactic halo is said to be self similar if its time evolution is such that the halo remains 
identical to itself except for an overall rescaling of its phase space density, and its size in 
spatial and velocity dimensions, by time dependent factors [19]. Under the assumption 
of self similar evolution with the added assumption of spherical symmetry, Fillmore and 
Goldreich [20] and Bertschinger [21] described the properties of galactic halos. This model 
was modified to include angular momentum by Sikivie, Tkachev, and Wang [2"2l 125] . For 
a recent review of the self similar model of the Milky Way halo, see [TH]. The self similar 
infall model predicts approximately flat rotation velocities far from the galactic center, in 
agreement with observations. The model is also consistent with the existence of a "core 
radius" observed in many galaxies. This is because most of the dark matter particles have 
angular momentum relative to the halo center and do not reach the central region, resulting 
in a depletion of dark matter particles relative to the spherically symmetric scenario [19]. 



A. Discrete flows and caustics. 

A prominent feature of the self similar model is the existence of cold discrete flows and 
dark matter caustics. Note that the existence of discrete flows and caustics is not a conse- 
quence of self similarity, but rather a consequence of Liouville's theorem. In [211 12S], it was 
argued that discrete flows and caustics should be a natural consequence of cold, collisionless 
matter. Each infall-outfall of dark matter produces an inner caustic and an outer caustic. 
Outer caustics are fold catastrophes that occur at the outer turnaround radii of particles, 
and appear as thin spherical shells surrounding galaxies. Their location is determined by 
the energy of the particles. Inner caustics occur near the inner turnaround radii of particles. 
Their location is determined by the magnitude of angular momentum, and their geometry 
is determined by the spatial distribution of the dark matter angular momentum field. For 
the special case of dark matter particles carrying a net rotation aligned with that of the 
baryons, the inner caustics are made up of elliptic umbilic catastrophes that resemble rings 
in the galactic plane [2H1 EZ]- Possible observational evidence for such ring caustics and 
for self-similarity of galaxies was found by [2H] who examined the rotation curves of spiral 
galaxies. The existence of caustics is relevant to the DAMA experiment since the velocity 
distribution in the vicinity of a dark matter inner caustic is dominated by the cold flow 
forming the caustic. In [29] , the effect of the dominant flow of the self similar model on the 
annual modulation signature was calculated. The effect of other cold streams on the recoil 
rate has been discussed by several authors (see for example [30H34"] ) . 

Table 1 (extracted from Table 1 of [H]) describes the first 40 flows in the self similar 
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infall model of the Milky Way halo. The first column is the fractional density contribution 
of each flow. The table is arranged in descending order of the flow density fraction. The 
dominant flow is assigned a fraction £. This flow is dominant because the associated inner 
caustic is close to the earth's location [THES] (In [H], a value of £ = 0.733 is adopted). We 
allow £ to be variable since the value of £ determines the peak of the annual modulation. 
£ is also very sensitive to the location of the closest inner caustic. The second and fourth 
columns give the maximum and minimum flow speeds relative to the earth. The third and 
fifth columns specify the time when the maximum and minimum occur respectively. 



B. Annual modulation. 

For a series of cold flows (i.e. ignoring the velocity dispersion, a valid assumption for 
WIMPs), the velocity distribution of WIMPs (relative to the halo) is: 

/flows (V) = $l£fl°ws,i<$(^- Vf,i), (13) 
i 

where £fl OW s,i represents the contribution of the i th flow to the local dark matter density. The 
mean inverse speed T(t, Q) can then be easily calculated for a flow vf. 



T(t,Q) = T^w-Mv^(t)\ - v min (Q)} 



l 



1 + — %(t) ■ V fe 
v fQ 



& b/0 - («©(*) •#/©)- Vmin(Q)] , (14) 



where Vf & and Vf Q are the flow velocities relative to the earth and sun respectively, and 9 
denotes the unit step function. %(i) is the velocity of the earth about the sun, and % is 
the velocity of the sun about the halo center. 

Fig. [2^a) shows the first 2 flows. Flow 1 is smallest in May and peaks in November. 
Flow 2 is smallest in July and peaks in January. Since Flow 1 is the dominant flow, it is 
instructive to look at this flow in detail. Let us choose a co-ordinate system in which the 
+x axis points towards the Galactic center, the +y axis points in the direction of Galactic 
rotation, and the +z axis points towards the north Galactic pole. In these co-ordinates, 
relative to the sun, Flow 1 has velocity [H] (but note that the co-ordinate system used in 
[T4] is different from ours): 

v hQ = 253.6 km/s (0.3549 x + 0.9345 y - 0.0276 z) (15) 

The velocities of the sun (about the halo center) and the earth (about the sun) in these 
co-ordinates are respectively (see for example [T2J |3Dl E2] , and references therein): 

v & = 233.3 km/s [0.0429 x + 0.9986 y + 0.0300 5] 
v @ (t) = 29.8 km/s [(0.9931 cos <p - 0.0670 sin 0) x 

+ ( 0.1 170 cos + 0.4927 sin <f>)y- (0.0103 cos + 0.8676 sin <p) z] , (16) 

where the angle <f>(t) = 2n (t — March 21)/365. We note that earth's velocity about the sun 
is most closely aligned with the sun's velocity about the halo center when <ft = 71°, which 
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The first 40 flows of the self similar infall model 
(from Table 1. of Q3|) 
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Aug 26 
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Mar 8 
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Sep 6 
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374 


Sep 10 


317 


Mar 12 


0.0550(1-4) 


311 


Oct 7 


257 


Apr 8 




0.0113(1-4) 


393 


Mar 25 


334 


Sep 24 


0.0550(1-4) 


355 


Jun 24 


325 


Dec 23 




0.0113(1-4) 


394 


Mar 24 


335 


Sep 22 


0.0550(1-4) 


382 


Dec 19 


322 


Jun 19 




0.0097(1-4) 


370 


Sep 6 


316 


Mar 8 


0.0324(1-4) 


365 


Mar 11 


307 


Sep 10 




0.0097(1-4) 


372 


Sep 7 


317 


Mar 9 


0.0243(1-4) 


380 


Mar 16 


321 


Sep 15 




0.0097(1-4) 


391 


Mar 27 


333 


Sep 26 


0.0227(1-4) 
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Jun 21 


400 


Dec 20 
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392 


Mar 26 


334 
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Mar 20 


337 


Sep 18 




0.0081(1-4) 


382 


Mar 29 


324 


Sep 27 


0.0146(1-0 


373 


Sep 13 


315 


Mar 15 




0.0081(1-4) 


383 


Mar 28 


325 


Sep 26 


0.0146(1-4) 


394 


Mar 21 


334 


Sep 20 




0.0065(1-0 


354 


Sep 2 


302 


Mar 4 


0.0129(1-0 


376 


Sep 11 


319 


Mar 13 




0.0065(1-0 


374 


Mar 31 


317 


Sep 29 


0.0129(1-0 


397 


Mar 23 


338 


Sep 21 




0.0049(1-4) 


635 


Jun 18 


579 


Dec 18 


0.0129(1-0 


529 


Jun 19 


477 


Dec 19 




0.0049(1-0 


644 


Dec 23 


598 


Jun 23 



TABLE I: The first 40 flows and their associated density fractions, from [14], in descending order 
of density contribution. The first column gives the fraction of the total density contributed by 
each flow. The dominant flow is assigned a density fraction 4- "max and " m in are the maximum 
and minimum flow speeds relative to the earth, seen at times i max and i m in- The recoil energy 
maximum observed by DAM A in the 2 — 6 keV ee range during May 17 < t < June 2 is obtained 
for 0.62 > 4 > 0.37. The mean DAMA best fit maximum on May 25 is obtained for £ = 0.47, 
while a peak on June 2 would correspond to a density fraction £ = 0.37. The flow densities 
published in |14j are obtained by setting 4 = 0.733 and multiplying by 231.8 x 10 - 26 gm/cm 3 . 
Note that the flow velocities in [H] are related to the velocities used here by the transformation 
x — > —x, y — > y, z — > —z. 

occurs around June 2. The two velocity vectors are most misaligned when = 251° which 
occurs six months later, around Nov 30. Flow 1 has speed relative to the earth: 

l*W*)l ~ v he ~ % [%(*) ■ «i,©] , (17) 

which is largest when = 225° (around Nov 5) and smallest when <fi = 45° (around May 7). 

Fig. [2]^b) shows v m i n for quenched recoils, for Na and I as a function of Qdet, for an 
assumed m x = 560 GeV (the value of m x is motivated by the data, as shown in the next 
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FIG. 2: (a) shows the first 2 flows of the self similar model, (b) shows v m i n for Na and I as 
a function of Qdet for a WIMP mass m x = 560 GeV, for quenched recoils. Also shown are the 
maximum and minimum velocities for Flow 1. (c) shows w m ; n for Na and I for quenched recoils, 
for different values of m x , for fixed Qdct = 6 keV ee - 
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section). Also shown are the minimum and maximum velocities for Flow 1, which occur on 
May 7 and Nov 5 respectively. For Q det < 7.4 keV ee , the flow is visible to the detector at all 
times of the year. We therefore expect a sinusoidal variation of dR/dQ with a peak in May 
(for this flow). For energies 7.4 keV ee < Qdet < 8.9 keV ee , the flow in invisible to Na during 
parts of the year. For energies 8.9 keV ee < Qdet < 10 keV ee , the flow is invisible to both 
Na and I during parts of the year, resulting in pronounced non-sinusoidal behavior. In the 
10 keV ee < Qdet < 12.1 keV ee range, the flow is completely invisible to Na and is seen by I 
only during some parts of the year near the flow maximum in November. For Qdet > 12.1 
keV ee , the flow is invisible to both Na and I, resulting in no recoils for this particular flow. 
Fig. [2|c) shows t> min as a function of WIMP mass, for Q det = 6 keV ee , for quenched recoils. 
For masses m x > 150 GeV, the flow is visible to both Na and I during all parts of the 
year. Thus for sinusoidal variation of dR/dQ at 6 keV ee , we must have m x > 150 GeV. We 
will see in the next section that more stringent bounds can be obtained. The recoil rate 
dR/dQ oc T(t, Q) oc 1/ \v\&\ provided the flow is visible to the detector at energy Q. 

It is interesting to contrast the self similar model with the Maxwellian. A Maxwellian 
halo is described by the distribution: 



fmaxi^Vwh) 



exp [-(v wh /v Q y 



7T" 



(18) 



where we have ignored the effect of the finite escape velocity. The subscript wh stands for 
"WIMP-halo" and indicates that the velocities are measured relative to the halo. Expressed 
relative to the detector, the (1-dimensional) speed distribution becomes: 



f(v) 



Vnv v eh 



( v + v eh Y 

\ v o J 



(19) 



implying a mean inverse speed 



-^max(^) Q) 



2Veh(t) 



erf ■ 



^min 

(Q)+Veh(t) 
V 



erf 



^min (Q) - v eh (t) 



(20) 



where v e h(t) = \v Q + v e (t)\. It is instructive to compute the angular dependence of the flux 
of dark matter particles on earth for the self similar infall model, and contrast it with the 
prediction of the isothermal halo, as done in [12] . Since the streams of the self similar model 
have a net velocity relative to the halo center, they do not all arrive from the direction of 
the sun's motion. Some streams arrive in directions above and below the galactic plane, 
but the densest streams (and in particular, the dominant, or big flow) are restricted to the 
galactic plane, in a direction nearly opposite to that of the sun's motion. In contrast, for the 
isothermal halo, the WIMP particles have no average velocity relative to the halo center, 
and the WIMP wind is due to the motion of the sun (and earth), implying a much larger 
flux in the direction of the sun's motion compared to the flux in the opposite direction. 

Also of note is the energy spectrum of recoils. In the self similar model, the velocity 
distribution is discrete. This means that at a fixed energy, a stream may or may not be 
visible, and the number of streams that contribute to the signal decreases as the recoil energy 
is increased. The mean inverse speed T given by Eq. [9] is therefore a series of steps, for 
the self similar model. The height of each step oc 1/v, where v is the speed of the flow 
under consideration (relative to the detector), and is therefore largest at times of the year 



9 



when the flow speed relative to the earth is the smallest. The edge of the step on the other 
hand oc v 2 , and is largest for the highest velocity flows. The modulation of the edge of 
each step is twice the modulation of the height, and has opposite phase. The net recoil 
rate integrated over energies depends not only on the flow speeds and densities, but also on 
the energy dependence of the cross section. As F 2 (Q) C 1 at high energies, the DAMA 
experiment is sensitive to recoils only for Qdct < 6 keV ee , and relatively low energy WIMP 
recoils contribute more than high energy recoils. 

For the self similar model, T peaks in May/ June (depending on the value of £) for 
large WIMP masses. For a fixed WIMP mass, T peaks in May/June at small energies, 
becoming non-sinusoidal, and possibly reversing phase at larger energies. In contrast, for 
the isothermal halo with a Maxwellian velocity distribution, T measured at a fixed recoil 
energy peaks in June for sufficiently small WIMP masses, but reverses phase and peaks in 
November for larger WIMP masses, as can be verified by expanding Eq. [20] in a Taylor series 
[39] . Conversely, for a fixed WIMP mass, T peaks in June at large energies, but reverses 
phase, peaking in November at very small energies. Given the very different nature of these 
two halo models, it is certainly surprising that the DAMA results may be fit to either model 
as we shall see in the following section (when the WIMP mass is not fixed by an independent 
measurement). 

IV. RESULTS 

We now compare the prediction of the self similar halo model with the DAMA results 
and obtain best fit values of mass and cross section. We set the dark matter density at our 
location p x = 0.3 GeV/cm 3 . We assume energy independent values for the quenching factor 
q.f. = 0.3 for Na and q.f. = 0.09 for I. We consider only the spin-independent cross section 
for WIMPs scattering elastically off a proton (or neutron) a p . The fraction of channeled 
recoils for Na and I as a function of energy has been calculated experimentally by [TS] • For 
a recent theoretical treatment, see [38] . Here we use the fit obtained by [36] to Figure 4 of 
|18j . for the channeling fractions for Na and I. 

We use Table 1 of [2] (summarized in Table 1 here) to obtain the flow velocities and 
densities. We use Table 2 published in [3] which gives the observed peak position in the 
recoil energy spectrum for different energies, in order to determine £. The amplitude of 
the recoil spectrum is presented in [1]. We fit the recoil rate at different energies to the 
observed amplitude of the annual modulation using Table 3 published in [36J which was 
extracted from Fig. 9 of [3]. We perform a minimum x 2 analysis using 36 energy bins and 
2 fitting parameters (m x and cr p ). The contours are obtained by plotting curves of constant 
X 2 = Xmin + ^x 2 - A% 2 is obtained by setting the area under the x 2 distribution equal to 
the required confidence level (C.L.), with n equal to the number of fitting parameters (see 



For n — 2, this simplifies to A% 2 = —2 log(l — C.L.), and we find A% 2 = 2.16, 5.99, and 
11.62 for 66%, 95%, and 99.7% confidence respectively. 



for e.g 



SS], I2Z!): 




(21) 
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A. Fitting the location of the maximum recoil rate. 

The maximum in the recoil spectrum measured by the DAMA experiment in the energy 
range 2-6 keV ee is t max = 144 (May 25) ±8 days. We vary £ to fit the DAMA phase. We 
find that £ = 0.47 fits t max = 144, implying that the dominant flow contributes 47% of the 
local dark matter density. Such a large contribution implies the existence of a nearby dark 
matter inner caustic (in [35], an observation of such a caustic is claimed). The observed 
maximum of 144 ± 8 days is obtained for 0.62 > £ > 0.37, with £ = 0.62 corresponding 
to a peak on May 17, while £ = 0.37 giving a peak on June 2. We note that the DAMA 
maximum picks Flow 1 as the dominant flow, while in [H] , the dominant flow is either Flow 
1 or Flow 2. £ is set to 0.47 for all our results. 



B. Best fit parameters. 

We fit our 2 free parameters m x and a p by minimizing 



36 



x 2 = E 



i=l 



^data,i ~~ ^4model,i(°"p) m \ 



(22) 



where the sum is over energy bins. Ad a ta,i is the measured amplitude for energy bin i and 
j 4modei,i( cr p ) m x) is the predicted amplitude for energy bin i for the assumed m x and a v . cr.j 
is the uncertainty in the measurement of Aj. We compute the amplitude v4 mode i as 

Anodcl = ^ I 'dQ 'det(£max) — dR/dQdct (^min)] , (23) 

in the energy region where A mo d e i is sinusoidal. For the self similar model, we use the DAMA 
best fit value of t max = 144. t m i n is set equal to 327. When a comparison with the isothermal 



halo is made, we use Eq. 18 with t> = 220 km/s, t max = 152, and t min = 335. 

Fig. [3] shows the expected recoil rate in the self similar infall model for 4 different 
energy bins. For energies Qdet < 6 keV ee , the recoil spectrum is qualitatively identical to 
(a). The expected sinusoidal variation with a maximum at t = 144 and a minimum at 
t = 327 is seen in (a). When the energy range is increased in (b), non-sinusoidal features 
start to appear since the flow velocity near t = 144 is less than the minimum velocity v m [ Q 
required to produce Na recoils at this energy (see Fig[2^b) for the dominant flow). As the 
energy is increased, t> m in increases, leading to the non-sinusoidal shapes seen in (c) and (d). 
The amplitude is negative in (d) indicating a phase reversal. DAMA does measure negative 
amplitudes at high energies, in particular the measurement in the energy range 9.5 — 10 keV ee 
is statistically significant. Smaller error bars are required before the negative amplitudes 
measured by DAMA can be treated as a physical effect. Negative amplitudes measured at 
high energies favor the self similar model (or a cold stream) and cannot be accommodated 
by the isothermal halo. The very small amplitudes seen in (c) and (d) due to the small value 
of F 2 (Q) make detection challenging. 

Fig. [4] shows the modulation amplitudes measured by the DAMA collaboration (open 
diamonds with error bars) at different energy bins. The prediction of the self similar model 
(m x = 570 GeV) is shown by the solid line, while the amplitudes expected for the Maxwellian 
models (m x = 78 GeV and m x = 12.4 GeV) are shown by broken lines. 
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FIG. 3: Recoil rate (mean subtracted) expected for the self similar infall halo model at different 
times of the year, for 4 different energy bins. The rate is a sinusoidal function of time when the 
flow velocity exceeds v m i n at all times of the year as in (a). For large energies [(c) and (d)], the 
flow is only visible to the detector during parts of the year, resulting in a non-sinusoidal pattern 
and phase reversal in (d). Ion channeling is included. 

Fig. [5] shows the 99.7% confidence contours in the mass-cross section plane for elastic spin- 
independent scattering, for the self similar infall model. The contours with solid lines include 
the effect of ion channeling while the contours with broken lines do not. The channeling 



effect is more important at lower WIMP masses. Using Eq. (22), the minimum value Xmm * s 



found to be 31.31/34 dof when the effect of channeling is included, and 34.93/34 dof without 
channeling at m x fa 570 GeV (dof stands for degrees of freedom). 

The contour for the self similar model is not closed. This is due to the fact that for 
WIMP masses m x much greater than the mass of Iodine, the minimum velocity 



/Qm N / Q 
Umin = V2^| "fc (24) 

becomes independent of WIMP mass m x . As a result, the recoil rate only depends on 
one parameter a p /m x , and not on a p and m x separately. Fitting the data to the single 
variable a p /m x , we find a best fit value of 0.059 ± 0.014 fb/TeV at the 95% level, with a 
Xmm °f 35.43/35 for the one parameter fit. The 66% contour (not shown) is closed between 
WIMP masses 330 GeV and 3.3 TeV. Shown for comparison are the contours for the simple 
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FIG. 4: Modulation amplitudes. Points are the DAMA/LIBRA measurements. The solid (red) line 
is the prediction of the self similar model, while the broken lines are drawn for the two Maxwellian 
models. 
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FIG. 5: Allowed regions in parameter space for elastic spin-independent scattering. The contours 
with solid lines include the effect of ion channeling, while the contours with broken lines do not. 
Shown are 3a results. The 3a contour for the self similar model is not closed at the high mass end 
because for m x much greater than the mass of an Iodine nucleus, w m ; n is nearly independent of 
m x . The la contour (not shown) is bounded between 330 GeV and 3.3 TeV. The allowed regions 
are ruled out by the CDMS and Xenon bounds (using data from 0). Also shown are the allowed 



regions for the Maxwellian halo of Eq. 18 Xmin ^ 

found to be 31.31/34 dof at m x = 570 GeV for 
the self similar model. For the Maxwellian we find Xmin = 30.6/34 dof at m x = 12 GeV, and Xmin 
= 26.41/34 dof at m x = 78 GeV. 



Maxwellian model of Eq. 18 The CDMS and XenonlOO exclusion contours are also shown 
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(CDMS and Xenon data from Fig. 5 of [6]). 



C. Adding a thermal component. 

Let us now modify our discrete sum over cold flows by adding a thermal component. This 
is similar to adding numerous flows that approximate a continuum to an experiment with 
finite energy resolution. We would naturally expect the innermost region of phase space to 
be unresolved to detectors while the outer region of phase space to be seen as a series of 
cold flows. We modify Eq. [13] as: 

Pf{v) — P [^flows/flows + (1 - £flows)/max] , (25) 

where £ now s is the contribution due of the 40 flows in Table 1, and / max is given by Eq. [18} 
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FIG. 6: Effect of adding a thermal component. The 3a contours with solid lines are for the 



distribution of Eq (25) with £ n ows = 0.7. The allowed region becomes larger as £fl OW s is reduced. 
The contours with broken lines are for different distributions (self similar or Maxwellian). Ion 
channeling is included. 

We choose £ now s = 0.7 so that the thermal component is comparable to the dominant flow. 
Fig. [6] shows the allowed 3cr contours. Also shown (in dashed lines) are the contours for 
the Maxwellian halo at low masses, and the self similar infall model for high masses (Note 



that the contours with solid lines are for the same f(v) (Eq. 25), and the contours with 
dashed lines are for different models). There exists a tiny region near m x = 12 GeV which 
is affected by some of the high velocity flows, but not by the dominant flow. The 3er contour 
is much smaller than for a pure Maxwellian halo, but becomes larger as £ n0 ws is decreased. 
Also cr p is about three times the value for a Maxwellian halo, as expected for £fl OW s = 0.7. 
Thus the cold streams of the self similar model do not help in bringing the DAMA result 
in agreement with other experiments. It is however important to note that the self similar 
model is consistent with a small WIMP mass provided a thermal component exists. The 
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other allowed region is at high WIMP masses m x > 250 GeV, as with the case of pure self 
similar infall, with slightly larger a p to compensate for the smaller value of £fl OW s = 0.7. The 
Xmin i s f° un d to be 30.64/34 dof near m x = 12 GeV and 30.05/34 dof near m x = 570 GeV. 
The x 2 function also has a minimum near m x = 48 GeV, but the large value (« 61) means 
it is disfavored at high significance. 

D. Comparison with the DAMA average plus unidentified background. 

The DAMA experiment measures only the modulation about the average recoil rate, not 
the average itself. This is because of the presence of a large background that does not 
modulate annually. Nevertheless, the DAMA collaboration has reported the sum of the 
average recoil rate and the unidentified background [4J. It is an important check that the 
average recoil rate predicted by a particular model is less than the average plus background 
reported by DAMA. Figure [7] shows the DAMA measurement (open diamonds) extracted 
from Fig. 1 of [1] compared with the self similar model (solid line) and the two Maxwellian 
models (broken lines). We see that the prediction of the self similar model is everywhere 
below the average plus background reported by DAMA, and is thus consistent with the 
observations. The Maxwellian model with m x = 78 GeV predicts an average rate that is 
too large in the first few energy bins. We will however not exclude that solution as our 
treatment of the isothermal halo is quite approximate. In [37], a more thorough analysis 
leads to conclusions similar to what we find here regarding the Maxwellian model with 
relatively large WIMP mass ~ 80 GeV. 




Energy (keV ee ) 

FIG. 7: Comparison with the measured average + background. The open diamonds (black) are 
the DAMA measurement of the average + background, extracted from Fig. 1 of [1]. The solid line 
(red) is the expected average recoil rate for the self similar model. We see that the line remains 
below the measured average + background at all energies. The Maxwellian models are shown by 
broken lines for m v = 12.4 GeV and m v = 78 GeV. 
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V. DISCUSSION. 



In this paper, we examined the DAMA annual modulation result in the context of the 
self similar infall halo model. We showed that the self similar model is in good agreement 
with the DAMA experiment for spin independent elastic scattering, with Xmm P er degree 
of freedom = 0.92(1.03) with(without) channeling, for WIMP masses exceeding 250 GeV at 
99.7% confidence. For large WIMP masses, the cross section-mass relation is approximately 
cr p /m x m 0.06(0.05) fb/TeV with(without) channeling. As in the case of the Maxwellian, 
the allowed region has been excluded by the CDMS and Xenon experiments. 

In Section II, we derived an expression for the expected recoil rate assuming a spin- 
independent cross section and elastic scattering. We then discussed the self similar infall 
model in Section III. We examined the speed of the dominant flow at different times of the 
year, and compared it to t> min for both Na and I (Fig. [2]). We then presented our results in 
Section IV. Fig. [3] shows the modulation amplitude (mean subtracted) expected for the self 

shows the best fit amplitudes compared 
shows the allowed regions in parameter 



similar model at four different energy bins. Fig 
to the DAMA/LIBRA measurement, while Fig. 
space with and without ion channeling taken into account. We then introduced a small 
thermal component and studied the effects (Fig. [6]). With a thermal component, there are 
two allowed regions in parameter space. The first is near m x = 12 GeV due to the thermal 
component. This region is only slightly affected by the flows. The other is near m x = 570 
GeV and is due to the cold flows. Finally, we verified (Fig. [7|) that the time averaged recoil 
rate predicted by the self similar model is consistent with the measurement of the average 
plus background reported by the DAMA collaboration. 

If the DAMA results are indeed correct, it is interesting to ask whether we can distinguish 
the different halo models. We have shown here that the self similar model has properties 
that distinguish it from the isothermal Maxwellian. Let us take a look at the differences: 



1. Negative amplitudes and measurements at low energies. 

As shown in Fig. [3j the recoil rate in the self similar infall model becomes non-sinusoidal 
at large energies, with the rate in November exceeding the rate in May/June. A convincing 
measurement of negative amplitudes (phase reversal) in the high energy bins would be 
consistent with the self similar infall model, but inconsistent with the isothermal Maxwellian 
halo. The DAMA result does include a statistically significant negative amplitude in the 9.5- 
10 keV ee bin. However most measurements at high energies are smaller than the uncertainty, 
and are thus not reliable. Moreover the scattering cross section is small at large energies, 
making the modulation difficult to measure. 

A phase reversal or non-sinusoidal behavior at low energies would be consistent with the 
Maxwellian, and inconsistent with the self similar model. As F(Q) ~ 1 at low energies, 
one may hope to use future data in these energy bins. The DAMA collaboration plans to 
add new photomultiplier tubes in order to achieve a lower energy threshold than the current 
threshold of 2 keV ee [3|. Let us suppose that future DAMA data includes measurements in 
the 1 — 1.5 keV ee and the 1.5 — 2 keV ee energy bins that are currently non-existent. Fig. [8] 
shows the modulation amplitudes in these energy bins. The Maxwellian model with m x = 78 
GeV is the only one that shows a negative amplitude (i.e. the recoil rate in June is less 
than the recoil rate in December at these energies). The Maxwellian model with m x = 12 A 
GeV shows a small positive amplitude, while the self similar model with m x = 570 GeV 
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shows the largest positive modulation amplitude. These models can thus be distinguished 
provided sufficiently small error bars are achieved. 
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FIG. 8: The modulation amplitude at energies below the current DAM A threshold of 2 keV ee 
(there is no data currently for Qdet < 2 keV ee ). The self similar model (solid lines) shows an 
increase in amplitude with decrease in measured energy, while the Maxwellian models (broken 
lines) predict a smaller amplitude with decreasing energy owing to phase reversal at very low 
energies. The Maxwellian with m x = 78 GeV shows a negative amplitude in the lowest energy bin 



2. Comparing the results of two different experiments. 

|4TH |4"T] have described a technique of calculating the WIMP mass using the results of 
experiments with two different target nuclei, by comparing the moments of the distribution 
function. The great advantage of this approach is that it does not assume a form for the 
distribution function, and can thus be called model independent. Once the WIMP mass 
is determined in this way, it is possible to place constraints on the distribution function. 
Similarly, a lower bound on the WIMP mass from accelerator experiments can also constrain 
the form of the velocity distribution in the solar neighborhood. The parameter space near 
m x = 80 GeV is extremely sensitive to the presence of streams. We have checked that when 
the streams of the self similar model contribute as little as 5% to the total dark matter 
density in the solar neighborhood (with the dominant stream contributing 2.35% of the 
total), the expected recoil rate disagrees with the DAMA result in the 2-2.5 keV ee bin at 
> 6a, for an assumed WIMP mass m x = 78 GeV. This is so because in this energy bin 
(and for this mass), the minimum velocity required to produce a recoil is comparable to the 
velocity of the dominant stream and as a result, the stream is visible only during certain 
months of the year. As a result, the modulation due to the stream rivals that of the entire 
halo at this energy (see for example, [30] for a discussion of this effect). Thus provided the 
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WIMP mass is determined to be ~ 80 GeV, the self similar model may be ruled out at high 
significance. Note also that only an experiment that is sensitive to the annual modulation 
will see this effect. The parameter space near ~ 10 GeV is far less sensitive to the presence 
of streams. This is because the minimum velocity required to produce recoils in the DAMA 
energy bins is often so large that all but the highest velocity streams are invisible. We find 
that stream fractions ^ 50% are admissible when m x ~ 10 GeV (see Fig. [6j. 

3. Directional sensitivity. 

For the Maxwellian halo, the WIMP particles do not have a definite direction, instead they 
have a large velocity dispersion. As a result, the dark matter particles come predominantly 
from the direction of the sun's motion [12]. This is not the case for the self similar infall 
halo. Since the flows of the self similar model are cold (i.e. non-thermal), the WIMP flux 
depends on the velocity vectors of the flows, and the most intense flows are in a direction 
nearly opposite to that of the sun's motion [TO] [12]. Thus a large WIMP wind due to 
the dominant flow can be easily distinguished from the prediction of the Maxwellian halo 
provided directional information is available. 

4- Measuring the average recoil rate in addition to the amplitude. 

As mentioned previously, the DAMA experiment does not measure the time averaged 
value of the recoil rate due to the large unidentified background. If this average can be 
measured, we have additional and complimentary information regarding the halo model. 
Fig. [9] shows the percentage modulation which is the ratio of the modulation amplitude to 
the average value of the recoil rate, for the three models. This ratio has the advantage that 
its energy dependence is entirely due to the halo model. The self similar model which has 
discrete streams produces a flat spectrum at low energies, since all streams are visible. The 
Maxwellian models on the other hand show an increase in the percentage modulation with 
energy. 

A detailed study of these ideas is left to future work. A major problem not considered 
in the present work is the incompatibility of the DAMA result with the null result of other 
experiments. As seen from Figs. [5]and[6j the self similar model does not help in bringing the 
DAMA result in agreement with other experiments. Inelastic scattering [I2HIE] has been 
suggested as a possible solution to this discrepancy. In the inelastic scattering scenario, the 
internal energy of the WIMP is altered, thus only a small fraction of the kinetic energy of 
the incoming particle is transferred to the target. Inelastic scattering prefers heavy targets 
and for a suitable mass splitting (i.e. energy difference between the lowest and excited states 
of the WIMP), one can observe recoils with Iodine, and very few (or none) with a lighter 
element such as Germanium. In this scenario, the CDMS bounds are considerably relaxed, 
while the bounds from Xenon, CRESST, and KIMS are still relevant. Further study of this 
scenario, as well as an analysis of the excess events claimed by the CoGeNT experiment 
[17] in low energy bins is left to future work. It is hoped that future dark matter detectors, 
particularly ones with directional sensitivity will be able to shed more light on the form of 
the local dark matter phase space distribution. 
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FIG. 9: Percentage modulation at different energies. The self similar model predicts a step like 
spectrum resembling a flat line at low energies. The Maxwellian models show a continuous variation 
with energy. 
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